function [unWarpPhase1] = GoldsteinUnwrap2D_r1(im_phase, im_mask)
%UNTITLED2 此处显示有关此函数的摘要
%   此处显示详细说明
im_mag   = abs(im_phase);                  %Magnitude image
mag_max = max(im_mag(:));
%indx1 = find(im_mag < 0.1*mag_max);  %Intensity = mag^2, so this = .04 threshold on the intensity
%im_mask(indx1) = 0;   % Don't mask at this point; wait until residues are computed
if(~exist('im_mask','var'))
    im_mask = ones(size(im_phase));          %Mask (if applicable)
end
% figure; imagesc(im_mag.*im_mask),   colormap(gray), axis square, axis off, title('GS Initial masked magnitude'); colorbar;
% figure; imagesc(im_phase.*im_mask), colormap(gray), axis square, axis off, title('GS Initial masked phase'); colorbar;

%% Compute the residues
residue_charge = PhaseResidues_r1(im_phase, im_mask); % Calculate phase residues (Does not use mask)
% figure; imagesc(residue_charge), colormap(gray), axis square, axis off, title('GS Phase residues (charged)'); colorbar;

%% Compute the branch cuts
max_box_radius=floor(length(residue_charge)/2);  % Maximum search box radius (pixels)
%max_box_radius=4  % Maximum search box radius (pixels)
if(~exist('max_box_radius','var'))
    max_box_radius=4;  % Maximum search box radius (pixels)
end
% BranchCuts() ignores residues with mask == 0, so keep the entire mask == 1
branch_cuts = BranchCuts_r1(residue_charge, max_box_radius, im_mask); % Place branch cuts
% figure; imagesc(branch_cuts),    colormap(gray), axis square, axis off, title('GS Branch cuts'); colorbar;

im_mask(branch_cuts) = 0;  % Now need to mask off branch cut points, in order to avoid an error in FloodFill
im_mag1 = im_mag.*im_mask; % Mask off magnitude == 0 points, so that they are not chosen for the starting point

%% Manually (default) or automatically identify starting seed point
if(1)  % Chose starting point interactively
    im_phase_quality = im_mag1;
    minp = im_phase_quality(2:end-1, 2:end-1); minp = min(minp(:));
    maxp = im_phase_quality(2:end-1, 2:end-1); maxp = max(maxp(:));
%     figure; imagesc(im_phase_quality,[minp maxp]), colormap(gray), colorbar, axis square, axis off; title('Phase quality map');
%     uiwait(msgbox('Select known true phase reference phase point. Black = high quality phase; white = low quality phase.','Phase reference point','modal'));
%  [xpoint,ypoint] = ginput(1);        %Select starting point for the guided floodfill algorithm
% fprintf("xpoint:%d  ypoind:%d \n",xpoint,ypoint);
xpoint = 19;
ypoint = 156;
    colref = round(xpoint);
    rowref = round(ypoint);
    close;                              %Close the figure;
else   % Chose starting point = max. intensity
    [r_dim, c_dim]=size(im_phase);
    im_mag1(1,:) = 0;                     %Set magnitude of border pixels to 0, so that they are not used for the reference
    im_mag1(r_dim,:) = 0;
    im_mag1(:,1) = 0;
    im_mag1(:,c_dim) = 0;
    [rowrefn,colrefn] = find(im_mag1 >= 0.99*mag_max);
    rowref = rowrefn(1);                  %Choose the 1st point for a reference (known good value)
    colref = colrefn(1);                  %Choose the 1st point for a reference (known good value)
end

%% Unwrap
if(exist('rowref','var'))
    im_unwrapped = FloodFill_r1(im_phase, im_mag, branch_cuts, im_mask, colref, rowref); % Flood fill phase unwrapping
else
    im_unwrapped = FloodFill_r1(im_phase, im_mag, branch_cuts, im_mask); % Flood fill phase unwrapping
end
unWarpPhase1 = im_unwrapped;
    
end

